/*
 * Copyright (c) 1997, 1998, 2003, 2006 Aleksandar Samardzic
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are
 * met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. The name of the author may not be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */

#include <assert.h>		/* for assert() function */
#include <stdlib.h>		/* for malloc() function */
#include "voxel_grid.h"

/* VoxelGrid_Init - VoxelGrid class constructor */
void
VoxelGrid_Init(this, pPrimitives, pSlabs)
     PVoxelGrid     *this;
     PList          *pPrimitives;
     PSlabs         *pSlabs;
{
	double          (*d)[2] = this->d;
	double          h[3],
	                l,
	                delta;
	Byte            slab,
	                power;
	DWord           numPrimitives = pPrimitives->numItems;
	PPrimitive     *pPrimitive;
	DWord           numCells,
	                cell;
	PListIterator   listIt;
	PBoundingVolume boundingVolume;

	for (slab = X; slab <= Z; slab++)
		d[slab][LO] = DBL_MAX, d[slab][HI] = -DBL_MAX;

	/* calculate bounding volume of scene */
	BoundingVolume_Init(&boundingVolume, pSlabs);
	for (ListIterator_Attach(&listIt, pPrimitives);
	     ListIterator_Next(&listIt, (void **) &pPrimitive);) {
		pPrimitive->CalcBoundingVolume(pPrimitive, &boundingVolume,
					       pSlabs);
		for (slab = X; slab <= Z; slab++) {
			if (boundingVolume.d[slab][LO] < d[slab][LO])
				d[slab][LO] = boundingVolume.d[slab][LO];
			if (boundingVolume.d[slab][HI] > d[slab][HI])
				d[slab][HI] = boundingVolume.d[slab][HI];
		}
	}
	BoundingVolume_Clean(&boundingVolume);

	/* number of grid cells should be roughly equal to number of
	 * primitives */
	power = 0, l = 1;
	for (slab = X; slab <= Z; slab++)
		if ((h[slab] = d[slab][HI] - d[slab][LO]) != 0)
			power++, l *= h[slab];

	/* calculate height of voxel grid cell */
	this->hCell = pow(l / numPrimitives, 1. / power);
	for (slab = X; slab <= Z; slab++)
		if (h[slab] == 0)
			h[slab] = this->hCell, d[slab][LO] -=
			    h[slab] / 2, d[slab][HI] += h[slab] / 2;

	/* calculate number of cells for each axis */
	for (slab = X; slab <= Z; slab++) {
		this->numCells[slab] = (DWord) ceil(h[slab] / this->hCell);
		h[slab] = this->numCells[slab] * this->hCell;
		delta = (h[slab] - (d[slab][HI] - d[slab][LO])) / 2;
		d[slab][LO] -= delta, d[slab][HI] += delta;
	}

	/* initialize voxel array */
	numCells =
	    this->numCells[X] * this->numCells[Y] * this->numCells[Z];
	Array_Init(&(this->cells), numCells, 0);
	for (cell = 0; cell < numCells; cell++)
		Array_AddItem(&(this->cells), NULL);

	/* voxelize all primitives in scene */
	for (ListIterator_Attach(&listIt, pPrimitives);
	     ListIterator_Next(&listIt, (void **) &pPrimitive);) {
		pPrimitive->rayID = 0;
		pPrimitive->Voxelize(pPrimitive, this);
	}
}

/* VoxelGrid_Clean - VoxelGrid class destructor */
void
VoxelGrid_Clean(this)
     PVoxelGrid     *this;
{
	PArrayIterator  arrayIt;
	PList          *pCellList;

	for (ArrayIterator_Attach(&arrayIt, &(this->cells));
	     ArrayIterator_Next(&arrayIt, (void **) &pCellList);)
		if (pCellList) {
			List_Clean(pCellList);
			free(pCellList);
		}
	Array_Clean(&(this->cells));
}

/* VoxelGrid_CellIntersect - find intersection of primitives from a voxel
 * grid cell with ray */
void
VoxelGrid_CellIntersect(this, cell, pParam)
     PVoxelGrid     *this;
     DWord          *cell;
     void           *pParam;
{
	struct {
		PRay           *pRay;
		PPrimitive     *pPrimitive;
		DWord           rayID;
	}              *pStruct = pParam;
	DWord          *numCells = this->numCells,
	    cellIndex;
	PList          *pCellList;
	PListIterator   listIt;
	PPrimitive     *pPrimitive;

	/* check if cell is in valid range */
	if (cell[X] < 0 || cell[X] >= numCells[X] || cell[Y] < 0
	    || cell[Y] >= numCells[Y] || cell[Z] < 0
	    && cell[Z] >= numCells[Z])
		return;

	/* calculate cell index */
	cellIndex =
	    cell[Z] * this->numCells[Y] * this->numCells[X] +
	    cell[Y] * this->numCells[X] + cell[X];

	/* return if cell is empty */
	if (!(pCellList = Array_GetItem(&(this->cells), cellIndex)))
		return;

	/* intersect ray with each primitive in cell */
	for (ListIterator_Attach(&listIt, pCellList);
	     ListIterator_Next(&listIt, (void **) &pPrimitive);)
		if (pPrimitive->rayID != pStruct->rayID
		    && pPrimitive->Intersect(pPrimitive,
					     (PLine *) pStruct->pRay))
			pPrimitive->rayID =
			    pStruct->rayID, pStruct->pPrimitive =
			    pPrimitive;
}

/* VoxelGrid_Intersect - find intersection of voxel grid with ray */
PPrimitive     *
VoxelGrid_Intersect(this, pRay)
     PVoxelGrid     *this;
     PRay           *pRay;
{
	static DWord    rayID = 1;
	struct {
		PRay           *pRay;
		PPrimitive     *pPrimitive;
		DWord           rayID;
	} params;

	params.pRay = pRay, params.pPrimitive = NULL, params.rayID =
	    rayID++;
	VoxelGrid_EnumerateCellsLine(this, pRay->p0, pRay->dp, &(pRay->t),
				     VoxelGrid_CellIntersect, &params);
	return params.pPrimitive;
}

/* VoxelGrid_Insert - insert item into grid cell */
void
VoxelGrid_Insert(this, cell, pItem)
     PVoxelGrid     *this;
     DWord          *cell;
     void           *pItem;
{
	DWord          *numCells = this->numCells,
	    cellIndex;
	PList          *pCellList;

	/* check if cell is in valid range */
	if (cell[X] < 0 || cell[X] >= numCells[X] || cell[Y] < 0
	    || cell[Y] >= numCells[Y] || cell[Z] < 0
	    && cell[Z] >= numCells[Z])
		return;

	/* calculate cell index */
	cellIndex =
	    cell[Z] * this->numCells[Y] * this->numCells[X] +
	    cell[Y] * this->numCells[X] + cell[X];

	/* get pointer to list of items in cell */
	pCellList = Array_GetItem(&(this->cells), cellIndex);
	if (!pCellList) {	/* there are no items in cell */
		pCellList = malloc(sizeof(PList));
		assert(pCellList != NULL);
		List_Init(pCellList);
		Array_SetItem(&(this->cells), cellIndex, pCellList);
	}

	/* insert item into cell */
	List_AddItem(pCellList, pItem);
}

/* VoxelGrid_EnumerateCellsLine - call function specified for each cell
 * along the line (according to [Aman87] */
void
VoxelGrid_EnumerateCellsLine(this, p0, dp, pt_max, pEnumFunction, pParam)
     PVoxelGrid     *this;
     double         *p0;
     double         *dp;
     double         *pt_max;
     void            pEnumFunction(PVoxelGrid *, DWord *, void *);
     void           *pParam;
{
	double          (*d)[2] = this->d;
	double          hCell = this->hCell,
	    _hCell = 1 / hCell;
	DWord           slab;
	PVector         S,
	                T;
	int             step[3];
	double          t_curr,
	                t_next[3],
	                t_delta[3];
	Byte            lo,
	                hi;
	double          t_grid[2],
	                t_slab[2];
	PVector         point;
	DWord           cell[3];
	Byte            drivingAxis;

	for (slab = X; slab <= Z; slab++) {
		S[slab] = p0[slab];
		T[slab] = 1 / dp[slab];
		step[slab] = Signum(dp[slab]);
		t_delta[slab] = step[slab] * this->hCell * T[slab];
	}

	if (p0[X] >= d[X][LO] && p0[X] <= d[X][HI] && p0[Y] >= d[Y][LO]
	    && p0[Y] <= d[Y][HI] && p0[Z] >= d[Z][LO] && p0[Z] <= d[Z][HI])
		t_curr = 0;
	else {
		t_grid[LO] = -DBL_MAX, t_grid[HI] = DBL_MAX;
		for (slab = X; slab <= Z; slab++)
			if (dp[slab] != 0) {
				t_slab[0] =
				    (d[slab][LO] - S[slab]) * T[slab];
				t_slab[1] =
				    (d[slab][HI] - S[slab]) * T[slab];
				if (t_slab[0] < t_slab[1])
					lo = 0, hi = 1;
				else
					lo = 1, hi = 0;
				if (t_slab[lo] > t_grid[LO])
					t_grid[LO] = t_slab[lo];
				if (t_slab[hi] < t_grid[HI])
					t_grid[HI] = t_slab[hi];
				if (t_grid[HI] < t_grid[LO])
					goto End;
			} else if (S[slab] - d[slab][LO] < 0
				   || S[slab] - d[slab][HI] > 0)
				goto End;

		t_curr = t_grid[LO];
	}

	Vector_LERP(point, p0, dp, t_curr);
	for (slab = X; slab <= Z; slab++) {
		cell[slab] =
		    (DWord) ((point[slab] - d[slab][LO]) * _hCell);
		if (cell[slab] == this->numCells[slab])
			cell[slab]--;
		t_next[slab] =
		    (dp[slab] !=
		     0) ? (d[slab][LO] + (cell[slab] +
					  ((step[slab] !=
					    -1) ? 1 : 0)) * this->hCell -
			   p0[slab]) * T[slab] : T[slab];
	}

	FOREVER {
		drivingAxis = (t_next[X] < t_next[Y]) ? X : Y;
		if (t_next[drivingAxis] == t_next[Z])
			t_next[drivingAxis] = t_next[Z];
		drivingAxis =
		    (t_next[drivingAxis] < t_next[Z]) ? drivingAxis : Z;

		(*pEnumFunction) (this, cell, pParam);
		if (t_next[drivingAxis] > *pt_max)
			break;

		cell[drivingAxis] += step[drivingAxis];
		t_next[drivingAxis] += t_delta[drivingAxis];
		if (cell[drivingAxis] < 0
		    || cell[drivingAxis] >= this->numCells[drivingAxis])
			break;
	}
      End:return;
}

/* VoxelGrid_EnumerateCellsCircleXY - call function specified for each
 * cell along the circle (according to [Andr97a]) */
void
VoxelGrid_EnumerateCellsCircleXY(this, center, radius, limit, pt,
				 pEnumFunction, pParam)
     PVoxelGrid     *this;
     double         *center;
     double          radius;
     double          limit;
     double         *pt;
     void            pEnumFunction(PVoxelGrid *, DWord *, void *);
     void           *pParam;
{
	Byte            slab;
	DWord           cell[3],
	                centerCell[2];
	double          x0,
	                y0;
	int             ri,
	                x,
	                y;
	double          rph,
	                rpx0,
	                rmx0,
	                rpy0,
	                rmy0;
	double          x0ph,
	                x0mh,
	                y0ph,
	                y0mh;
	double          d,
	                d0;

	for (slab = X; slab < Z; slab++)
		centerCell[slab] = (int) (center[slab] - 0.5);
	x0 = (center[X] - 0.5) - centerCell[X];
	y0 = (center[Y] - 0.5) - centerCell[Y];
	cell[Z] = (int) center[Z];

	ri = (int) (radius + 0.5);

	if (ri <= 1) {
		/* circle with small radius */
		d0 = 2 * radius + 1;
		for (x = -2; x <= 2; x++)
			for (y = -2; y <= 2; y++) {
				d = (x - x0) * (x - x0) + (y - y0) * (y -
								      y0) -
				    radius * radius;
				if (d >= 0 && d < d0)
					if (d < limit)
						cell[X] =
						    centerCell[X] + x,
						    cell[Y] =
						    centerCell[Y] + y,
						    (*pEnumFunction) (this,
								      cell,
								      pParam);
			}
	} else {
		/* general case */
		rph = radius + 0.5;
		rpx0 = radius + x0, rmx0 = radius - x0, rpy0 =
		    radius + y0, rmy0 = radius - y0;
		x0ph = x0 + 0.5, x0mh = x0 - 0.5, y0ph = y0 + 0.5, y0mh =
		    y0 - 0.5;

		/* starting point computation */
		y = 0, x = ri;
		if ((d =
		     ((x - x0) * (x - x0) + y0 * y0 -
		      radius * radius) / 2) < (radius - (x - x0))) {
			d += (x++) - x0mh;
			if (d < limit)
				cell[X] = centerCell[X] + x, cell[Y] =
				    centerCell[Y] + y,
				    (*pEnumFunction) (this, cell, pParam);
		} else if (d >= 0 && d < rph) {
			if (d < limit)
				cell[X] = centerCell[X] + x, cell[Y] =
				    centerCell[Y] + y,
				    (*pEnumFunction) (this, cell, pParam);
		} else {
			d += x0ph - (x--);
			if (d < limit)
				cell[X] = centerCell[X] + x, cell[Y] =
				    centerCell[Y] + y,
				    (*pEnumFunction) (this, cell, pParam);
		}

		/* generation of quadrant 1 */
		while (x > 0)
			if (d < rpy0 - y) {
				d += (y++) - y0mh;
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else if ((d += x0ph - (x--)) >= 0) {
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else {
				d += (y++) - y0mh;
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			}

		/* boundary points between quadrant 1 and quadrant 2 */
		if (d < rpy0 - y) {
			d += (y++) - y0mh;
			if (d < limit)
				cell[X] = centerCell[X] + x, cell[Y] =
				    centerCell[Y] + y,
				    (*pEnumFunction) (this, cell, pParam);
		}

		/* generation of quadrant 2 */
		while (y > 0)
			if (d < rmx0 + x) {
				d += x0ph - (x--);
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else if ((d += y0ph - (y--)) >= 0) {
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else {
				d += x0ph - (x--);
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			}

		/* boundary points between quadrant 2 and quadrant 3 */
		if (d < rmx0 + x) {
			d += x0ph - (x--);
			if (d < limit)
				cell[X] = centerCell[X] + x, cell[Y] =
				    centerCell[Y] + y,
				    (*pEnumFunction) (this, cell, pParam);
		}

		/* generation of quadrant 3 */
		while (x < 0)
			if (d < rmy0 + y) {
				d += y0ph - (y--);
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else if ((d += (x++) - x0mh) >= 0) {
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else {
				d += y0ph - (y--);
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			}

		/* boundary points between quadrant 3 and quadrant 4 */
		if (d < rmy0 + y) {
			d += y0ph - (y--);
			if (d < limit)
				cell[X] = centerCell[X] + x, cell[Y] =
				    centerCell[Y] + y,
				    (*pEnumFunction) (this, cell, pParam);
		}

		/* generation of quadrant 4 */
		while (y < 0)
			if (d < rpx0 - x) {
				d += (x++) - x0mh;
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else if ((d += (y++) - y0mh) >= 0 && y) {
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			} else {
				d += (x++) - x0mh;
				if (d < limit)
					cell[X] =
					    centerCell[X] + x, cell[Y] =
					    centerCell[Y] + y,
					    (*pEnumFunction) (this, cell,
							      pParam);
			}
	}
}
